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Abstract 

Numerical stochastic integration is a powerful tool for the investigation of quan- 
tum dynamics in interacting many body systems. As with all numerical integration 
of differential equations, the initial conditions of the system being investigated must 
be specified. With application to quantum optics in mind, we show how various com- 
monly considered quantum states can be numerically simulated by the use of widely 
available Gaussian and uniform random number generators. We note that the same 
methods can also be applied to computational studies of Bose-Einstein condensates, 
and give some examples of how this can be done. 

PACS numbers: 02.60.Cb, 02.60.Jh, 42.50.-p 
Key words: Numerical simulation, quantum states, quantum dynamics. 



1 Introduction 



The theoretical study of non-equilibrium quantum many-body dynamics is a 
growing area, especially since the experimental achievement of trapped Bose- 
Einstein condensates. Many of the methods used for theoretically investigating 
condensates have been adapted from theoretical quantum optics pQ, with vary- 
ing degrees of success. One particular approximation technique that proved ex- 
tremely successful in quantum optics is hnearisation of the fluctuations about 
solutions of the classical equations of motion. This technique, if used appro- 
priately, is a very powerful tool for the calculation of the steady-state spectra 
of intracavity parametric processes p]. However, in a dynamically evolving 
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system, or one operating near phase transitions or critical points, this method 
can give incorrect answers [31H] . The vahdity of the approximation depends on 
three conditions. The first of these is that the solution of the classical equa- 
tions is the same as the mean-field solution of the full quantum equations. The 
second and third are that the fluctuations about these solutions are in some 
sense small and that they can be represented as Gaussian, so that moments 
of higher than second order vanish. In the study of trapped Bose-Einstein 
condensates, the Hartree-Fock-Bogoliubov (HFB) method is a closely related 
approximation [3], and therefore needs to be used with the same care as the 
linearised fluctuation approximation in quantum optics. 

When these conditions are not met, there are still a number of ways to pro- 
ceed. In some very rare cases it may be possible to solve directly either a 
master equation for the density matrix, or even the Heisenberg equations of 
motion for the actual system operators. However, the most interesting quan- 
tum dynamics are not generally restricted to such cases. One set of methods 
which has been very successful is the phase-space representations originally 
used to develop stochastic differential equations in quantum optics ^. These 
allow common classes of quantum Hamiltonians to be mapped via master and 
Fokker-Planck equations onto stochastic differential equations. In some cases 
the Fokker-Planck equation may be solved directly for a pseudoprobability dis- 
tribution which then allows for the calculation of operator moments [T|7] . Once 
again, these cases are rare and can often only be solved in the steady-state 
regime. The method of choice if we wish to obtain dynamical quantum infor- 
mation is then to numerically integrate the stochastic equations of motion. As 
with any numerical analysis of differential equations, this then requires that 
the initial conditions be specified, as these can have marked effects on the sub- 
sequent dynamics, in both optical p[9] and interacting atomic and molecular 
systems |10|ll|12|13lll4lll5|16j . In what follows we will begin with a brief out- 
line of the theory behind the phase-space representations and then show how 
to numerically simulate some of the more common and useful initial quantum 
states of optics and condensed atom physics, in both the positive-P [T7] and 
Wigner representations [T8] . 



2 Phase-space representations of the density matrix 

Phase space techniques are a powerful tool to investigate the full quantum 
dynamics of interacting quantum systems in cases where it is impractical to 
solve either the Heisenberg equations of motion or the master (von Neumann) 
equation for the density matrix. Instead of working with operators or den- 
sity matrices, they allow us to work directly with classical c-number variables, 
which are amenable to manipulation on available computers. Perhaps more 
importantly, the complexity of the computation scales with the number of in- 
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teracting modes rather than with the size of the Hilbert space, which is often 
completely intractable. In fact, a single-mode quantum calculation has been 
performed using these methods for the order of 10^'^ interacting quanta |19j , 
which would be completely out of the question using other methods. There 
are a number of phase-space representations, among them being the Wigner 
representation [18], the Glauber-Sudarshan P representation [20121] . the Q 
representation [22] (sometimes known as the Husimi representation), the com- 
plex P representation [17J, and the R representation [20]. The most useful for 
numerical work are the positive-P and truncated Wigner representations [23], 
the latter being an approximation to the full Wigner representation. 



2.1 Truncated Wigner equations 



Historically, the first of these phase-space representations was the Wigner 
representation [18] , which was formulated as a pseudoprobability function for 
the position and momentum of a particle. Mathematically, the quadrature 
phase amplitudes of quantum optics are completely equivalent to position and 
momentum, so that the Wigner function is a frequently used tool for describing 
nonclassical states of bosonic fields. Quantum mechanical expectation values 
for operator products expressed in symmetrical order are found naturally in 
the Wigner representation as classical averages of the corresponding Wigner 
variables. As an example, making the correspondence between the single-mode 
annihilation operator a and the complex Wigner variable a, we find that 

-^=\(^ci^a + aa^) = N + ]^, (1) 

where N is the number of quanta in the mode. Given a general Hamiltonian 
which is some combination of bosonic creation and annihilation operators, if, 
we find the von Neumann equation as 

^^f=[^,P], (2) 

from which the equation of motion for the Wigner function, W , is found using 
the correspondence rules. 



Following the standard methods ^24j, as long as the equation found by the 
above procedure has no derivatives of higher than second order, it can be 
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mapped onto a set of stochastic differential equations for tlie variables a and 
a*. Unfortunately, all interesting problems result in derivatives of third order 
or more and, although methods exist for mapping the resulting generalised 
Fokker-Planck equations onto stochastic difference equations ^23], these are 
not very useful in practice. A common practice is to truncate the partial 
differential equation for the Wigner function at second order, often justified by 
claiming that the effect of these terms is small. This procedure may be formally 
justified by requiring the system modes to be highly occupied, and results in 
stochastic differential equations in what is known as the truncated Wigner 
representation. If there are no second order derivatives, the resulting equations 
are regular and quantum noise enters via the initial Wigner distribution for 
the variables. In optical problems, this then becomes functionally equivalent to 
stochastic electrodynamics [26] and has been shown to give misleading results 
in some cases [27P5] . This approximate method has also been used with some 
success in the study of Bose-Einstein condensates [2^30113 1(52] and is closely 
related to "classical field methods", including the stochastic Gross-Pitaevski 
equation [33|34|35ll36] . The appropriate initial states to use in the truncated 
Wigner equations are exactly the same as those that would be used in a full 
Wigner representation, with the approximations entering into the equations 
of motion. 

2.2 Positive-P representation 

The Glauber-Sudarshan P representation [2n][2T] is another representation of 
the density matrix in terms of coherent states and gives averages of the phase- 
space variables which are equivalent to normally-ordered operator expectation 
values. 



As photodetectors naturally measure normally-ordered averages, this would at 
first glance seem to be an extremely useful representation. It does, however, 
have two serious drawbacks. The first is that it is difficult to represent any 
state which is "more quantum" than a coherent state, as these do not possess 
positive and analytic P-functions. Although a P-function can be written for 
any quantum state in terms of generalised functions [37J, it is difficult to see 
how to sample these numerically. The second drawback arises when we con- 
sider the P-representation Fokker-Planck equation, found using the operator 
corresp ondences 




(4) 




(5) 



4 



It is readily seen that, for any interesting problem, the resulting Fokker-Planck 
equation will not have a positive-definite diffusion matrix and therefore will 
not be able to be mapped onto stochastic differential equations. The positive- 
P representation [.17j was developed to circumvent this problem by using a 
doubled phase space. For Hamiltonians which lead to derivatives of no higher 
than second order, this results in a Fokker-Planck equation which always has 
a positive-definite diffusion matrix and therefore can always be mapped onto 
stochastic differential equations. The price which has to be paid is that, instead 
of having a and a* as complex conjugate variables, the variables correspond- 
ing to this pair become independent. These are written in various ways, but in 
this article we will write the pair as a and a"*", and the appropriate equations 
can be found by naively using the P representation correspondences of Eq. [5] 
and then substituting a"*" for a* . The independence of the variables can cause 
serious stability problems with the numerical integration, but for problems 
where the integration converges, the positive-P representation is an extremely 
powerful theoretical tool [38] • As a final remark, we note that a method has 
been developed for mapping Hamiltonians which would give higher than sec- 
ond order derivatives in a generalised Fokker-Planck equation onto stochastic 
difference equations [39j, which is useful for analysing processes such as third 
harmonic generation [ID] and others which go beyond the common three and 
four-wave mixing processes of quantum optics and trapped ultra-cold gases. 



3 Representation of quantum states 

As always with the numerical integration of differential equations, initial con- 
ditions must be specified. In this section we will describe how this can be 
done for a number of quantum states which typically arise in quantum op- 
tics and BEC problems. We note here that we do not specify the type of 
stochastic differential equation, and all our results may be used as initial con- 
ditions for either regular (common in the truncated Wigner for systems with 
quadratic Hamiltonians) or stochastic (common in positive-P and truncated 
Wigner for damped systems) differential equations. For simplicity, we will con- 
sider single-mode representations of the states, with a multi-mode extension 
being relatively straightforward. We will also show how our choices for these 
states represent the correct quantum statistics, in terms of both quadrature 
and intensity variances. Note that we use quadrature definitions such that 





(6) 



so that the coherent state variances are equal to 1. 
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3.1 Coherent states 



The simplest initial condition to model in the positive-P representation is a 
coherent state, |/?), defined by a|/9) = When the Glauber P-function for 
a given state is well-behaved, we may also use this as the positive-P function. 
This means that a coherent state can be represented by the pseudoprobability 
distribution P(a,a"'") = 5{a — [3)5{a^ — /?*). Numerically, this means using 
the same complex conjugate pair, [3 and (3* , for each of the stochastic trajec- 
tories. This state is an appropriate choice to represent a well-stabilised laser 
output and is also commonly used as an initial condition for single-mode BEG 
analyses. 

The Wigner function for the coherent state is a Gaussian, 

= -exp f-2|a . (7) 

This is simply represented numerically by choosing the initial condition for 
each trajectory as 

1 

a = /5 + 2 (^1 + ^^2) , (8) 
where the rjj are sampled from a real normal Gaussian distribution, such as 
is given by the Matlab function randn. These have the correlations 77J = 
and rjjfjk = Sjk, where the overline denotes an average over many samples. 
It is readily shown that a = (3 and |a|2 = + | and that the quadrature 
variances both give 1, as required. 

3.2 Thermal or chaotic states 

These states, which can be used to represent, for example, a mechanical oscil- 
lator in a thermal bath [41j, have a particularly simple P-function, with 

m = 4exp(-|/5|Vrr), (9) 
irn 

where n is the average number present in the mode pj. We see immediately 
that, as expected, there is no phase information in this state. The appropriate 
distribution can be sampled from a normal Gaussian distribution multiplied 
by ^/w and a random phase term, 

a = \fTiri X exp {2t:%Q , (10) 

where 77 is again a normal Gaussian variable and C is uniformally distributed 
on [0,1]. Gaussian and uniform variables may be easily sampled using, for 
example, the Matlab functions randn and rand respectively. Numerical checks 
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of this distribution show that, with sufficient samples (> 10^), it reproduces 
well the intensity, n, and the quadrature variances, V{X) = V{Y) = 2n + 1. If 
the state is a mixture of coherent and chaotic states, i.e. a chaotic state with 
a coherent displacement, the P-f unction is written as 

P{P) = ^exp{-\(3 - (3o\Vn), (11) 

where f3o is the coherent displacement. This may be easily sampled as 

a = Pq + V^t] X exp (27riC) , (12) 
where the random variables are as in Eq. [TOl 

The Wigner function for the chaotic state may be found as a convolution of 
the P-function with a Gaussian of standard deviation one-half. In this case, 
where the P-function is itself a Gaussian, this results in a broader Gaussian, 
which can be sampled via 

a = ,Jn + l/2r] X exp (27rO , (13) 

where the random terms are the same as in Eq. [101 Samples from this distribu- 
tion reproduce well the required intensities and quadrature variances, always 
remembering that the Wigner distribution represents symmetrically ordered 
operator products. 



3.3 Squeezed states 

Squeezed states are those in which the variance of one quadrature of the field 
is below the coherent state value. Due to the Heisenberg uncertainty principle, 
which requires that V{X)V{Y) > 1, this means that the variance in the con- 
jugate quadrature must be greater than that of a coherent state [12]. They can 
now be routinely produced in the laboratory and are a useful resource for such 
things as nonclassical pumping of parametric processes, outcouphng squeezed 
atom laser beams [13] and quantum-limited measurement [Hj. Theoretically, 
a minimum uncertainty squeezed state with V{X)V{Y) = 1 is defined by the 
action of the squeezing operator, 

^(e)=exp[(e*)V/2-e2(at)V2], (14) 

on the vacuum, followed by the coherent displacement operator, 

D(r/)^(e)|0) = |r/,e), (15) 
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where e = re^"'^ |T]. The role of the squeeze factor, r, is apparent when we look 
at the quadrature variances (setting = for convenience), 

V{X) = e-2^ V{Y) = e^^ (16) 

Another point to note is that the squeezing process adds quanta to the mode, 
so that 

(r], e|d^a|r7, e) = |?7p + sinh^ r. (17) 

We will demonstrate here how to develop a numerical simulation of squeezed 
states using a canonical expression for an arbitrary positive-P function [17j. 
Given a density matrix p, a particular form of the positive-P function is 

P(a,a+) = ^|((« + (a+)*)/2|p|(a + (a+)*)/2)|Vl"-("")*l^/^ . (18) 

We now use the linear transformation 

p=(a + (a+)*)/2, 7 = (a-(«+)*)/2, (19) 
which has the Jacobian 



(20) 



a(p^,Py,7^,7y) 

so that the normalised distribution in terms of the new variables becomes 

P„.,).M£:^. pi) 

TT TT 

The displacement property 

D\ii)D{r]) = D\fx - 7]) exp [-ilm{fxri*)], (22) 
can then be used to obtain 

(/x|r/,e) = (0|Dt(^_^)5(,)|o)expHlm(pr/*)]. (23) 
It is now convenient to make another change of variables by setting 

fi-r] = e'^v, (24) 
and also make use of the normally ordered form of the squeeze operator, 

S{r,(j)) = (cosh r)""*^/^ exp (j"^'^^'^^ ^-^P [— ln(coshr)a^a] exp (j^^^^^ (^5) 
where F = e^*'^ tanh (r) . Making these substitutions in Eq. [23], we find 

K"!"-^*!'" ^STw • 



8 



and finally arrive at the separable Gaussian form 

P(i/, 7) = , , (27) 

V vre cosh r v ne^ cosh r vr 

which can be sampled using standard methods. We sample Eq. [27] and invert 
to find the appropriate random variables for the positive-P distribution initial 
condition as 



a = e''t'u + -f + r], (28) 
a+ = e-'^u* -'y* + 7]*, (29) 

where, given Gaussian random variables satisfying rTJ = and njuj: = Sjk, 



7 = -^(ni + m2), (30) 

/e~''coshr . /e'^coshr 
iy = d ns + td n4. (31) 

This distribution can be checked numerically and gives the appropriate values 
for the quadrature variances and intensities, with more samples needed for 
accuracy as the squeezing parameter becomes larger. 

The squeezed states are very simply modelled in the Wigner representation, 
by deforming the coherent state distribution in the appropriate manner. For 
a squeezed state with coherent dispalcement rj, this is done by sampling 

« = ^ + ^(Cie-^ + <2e'-), (32) 

with the ( being real normal Gaussian variables. A numerical check of this dis- 
tribution shows that it also reproduces the analytically calculated intensities 
and quadrature variances very accurately. 



3.4 Number or Fock states 



The Fock state is a quantum state with a fixed number of quanta. In optics, 
for example, the decay of an excited two-level atom can give a one photon 
Fock state. Here we will give a method which was previously used to model 
spontaneous emission from single-atom bosonic states [15] in the positive-P 
representation. The Fock state with n quanta has density operator p = \n){n\, 
which we will again sample using Eq. [TSl Introducing the new variables 

/i = and 7 = , (33) 
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we find the separable expression 



e-l^l' |;u|2"e-l^l' 
e-l^l'r(|/ip,n + l) 

TT TT 



(34) 



where 



r(^'^)= (35) 



is the Gamma distribution. Once again 7 = (ni+in2)/ is easily sampled via 
standard methods, while the Gamma distribution may be easily and efficiently 
sampled using a method given by Marsaglia and Tsang [16] to give z = 
and /i = ^/z e*^, where 6 is uniform on [0, 2tt). We then invert to find 

a = /i + 7 and = n* — 7*, (36) 

which are now correctly distributed to represent the Fock state \n). Numerical 
checks once again show that the intensities and variances are represented well 
if enough samples are taken. As seen in Ref. [15], the use of this method for 
the initial state also leads to analytically known dynamics being reproduced. 

To model a Fock state in the Wigner representation, we adapt an approxi- 
mation developed by Gardiner et al. [33] which allows us to approximately 
represent these without having to deal with negative pseudoprobabilities. The 
Wigner function for the Fock state |A^) is 

W{a, a*) = 2^ exp(-2|anL^(4|an, (37) 

TT 

where Lat is the Laguerre polynomial of order N. This distribution is oscil- 
latory and can obviously be either positive or negative, so cannot be easily 
simulated numerically. However, in the large N regime Gardiner has made 
the observation [33] that the cumulative distribution behaves very like a step 
function centered at |ap = N. This distribution can then be approximated 
by a Gaussian which gives the right moments for the mean and variance and 
approximates the higher moments well. The appropriate distribution is 

-"(-^)^\/!-(-^^^^|^)' 

where we have taken a = A/ne*^, with 6 uniform on [0,27r). The ffist three 
moments of this distribution are 
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a*a = N + - (39) 
^^=(^iV+iy + l (40) 

3 



IV 3 1 



a*3a3=(^iV+_j +_^N+-y (41) 

so that mean and variance are in agreement with what we expect analytically. 
We can now show that such an approximation in fact generates all moments 
of (!37|l . up to a correction of order 1/iV^, which is negligible for large N . 

Using the differential recursion relation for the Laguerre polynomials, 

x^Ln{x) = N{Ln{x) - Ln-i{x)), (42) 
ax 

a recursion relation for arbitrary even moments can be found. Writing 

(a*-a™)^ = 2^ ^— / d^a exp (-2|a|2)L;v(4|a|') (43) 

we can use to find 



, , N + m , , 

(q;*"^Q!"^)jv = — ^ — (a*™-ia™-i)jv + y (Q!™-^Q!'"~^)iv-i- (44) 
The first three moments are 

(tf^)^ = Ar+ 1/2, (45) 
(^^)^ = (AT + 1/2)2 ^ (4g) 

(^^3^)^ = (AT + 1/2)3 + ^(Ar + 1/2), (47) 

and comparison with (1391) shows that the Gaussian approximation is exact for 
the mean and variance, and accurate to 0(1/A^2) for m = 3. It is easily shown 
by induction on m that the exact moments satisfy 

{a*'^a^)^ = (A^ + 1/2)"^ + 0(Ar"-2), ^^g^ 

so that the correction is always of order l/A^^ relative to the leading term. 
Returning to the Gaussian approximation, we see that it gives 
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/ 2 

J- / dz{z + N + 1/2)™ exp (-2^2) 

V TT J-oo 



' 2 /"oo 

dz 

TT J-oo 



(iV + 1 /2)'" + m(Ar + 1 /2)"'~^z + OiN""-') exp {-2z') 



= {N +1/2)"" + 0{N"'~'), (49) 

which will be an adequate description of the number state statistics for large 
N. 

To simulate this distribution numerically, consider the choice 

a = p + qr], (50) 

where 77 is a normal Gaussian random variable, and p and q are yet to be 
determined. As we are using a Gaussian approximation, it is sufficient to 
reproduce the first two moments of a^. (Note that a is a real variable here, with 
the phase distribution to be added later). We need to reproduce = 1/2 
and ^ = (A^ + 1/2)^ + 1/4. The choices 

p = ^(2N+l + 2 ViV2 + iv) , (51) 

and 

reproduce the required distribution to a high degree of accuracy. The a thus 
chosen are then multiplied by the factor exp(2i7r^), where ^ is randomly chosen 
from the uniform distribution [0, 1). 



3.5 Crescent states 



The crescent state is given this name because its Wigner contours are sheared 
in phase-space due to a x^^^ (Kerr) nonlinearity, and is consistent with quan- 
tum states which have been proposed for trapped Bose-Einstein condensates ji7|48ll49j . 
where s-wave scattering is equivalent to a Kerr nonlinearity. In this section we 
will treat the Wigner distribution first, because we will use the corresponding 
Q-distribution to sample a positive-P distribution for the crescent state. 

In the Wigner representation a sheared state with coherent displacement ao 
is simulated by beginning with the squeezed state representation given pre- 
viously, and transforming this by a factor exp(ig?73), where q is the shearing 
factor, 



a 



"0 + ^ ('7ie + ^?7e' 



e''^''\ (53) 
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The real noise terms have the correlations 



Vj = 0, r]ir]j = 6ij. 



(54) 



Numerical checks of distributions produced using these methods again show 
that they give the expected values for average numbers and quadrature vari- 
ances. 



The Q-distribution Q{fi,fi*) can be simulated as a simple broadening of the 
Wigner distribution, so that we sample 



(55) 



We can then make use of Eq. f[T8|) to construct samples of the corresponding 
positive-P distribution: transforming to the variables 



/i=(a + (a+)*)/2, 



7 = («-(«+)*)/2, 



(56) 



we have 



. . (57) 

SO that given the Q-function samples (1551) . and 7 = (ni + in2)/V^, we have 
the crescent state sampling for the positive-P distribution 



a = /i + 7 and 



/i -7 



(5J 



3.6 Interacting many body states 



Generating the many body states appropriate for ultra-cold Bose gas simu- 
lations is difficult to illustrate within the single mode approach taken in this 
article. We note, however, that a number of Wigner states which are commonly 
used in modelling Bose-Einstein condensates may be sampled relatively easily 
and we provide a brief outline of currently available methods. We refer the 
reader to Ref. [36] for a detailed review of Wigner sampling methods for Bose 
gases. The methods may be summarized as follows 

1) Coherent state. — At zero temperature, a first approximation to the state 
of a BEG is a coherent state. In a continuous field theory the modes 
orthogonal to the condensate must necessarily contain vacuum noise in 
the Wigner representation. A simple and effective means to construct a 
coherent state is to simply add this noise, corresponding to half a quanta 
per mode, to the appropriate mean field solution of the GPE. 

2) Bogoliuhov state. — Steel et al. [29] demonstrated that an improved ap- 
proximation at low temperatures, the Bogoliubov state, may be readily 
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constructed from a stationary solution of the Gross-Pitaevskii equation 
once the Bogohubov modes are known. This approach has also been ex- 
tended to include U{1) symmetry constraints imposed by number conser- 
vation by Sinatra et al. [5U|31f32j . 

3) Adiabatic mapping. — Polkovnikov and Wang [35j used the quantum adi- 
abatic theorem to show that interacting states at zero temperature may 
be obtained by sampling the appropriate noninteracting state and then 
adiabatically ramping up interactions to the desired final value. 

4) High temperature states. — In the vicinity of a first approximation for 
the Wigner distribution is given by the non-interacting Bose-Einstein 
distribution. This may be used as a starting point for evolution according 
to the stochastic Gross-Pitaevskii equation [Ti|15f33IIM] which evolves 
the Bose field toward a sample from the grand canonical ensemble for the 
many body system. 



4 Conclusions 

In conclusion, we have shown how to take numerical samples in phase space of 
some of the most commonly appearing quantum states in quantum and atom 
optics. These techniques are important when we wish to investigate the effects 
of different initial states on dynamical quantum processes and possess straight- 
forward generalisations to many-mode problems. These methods will become 
more useful and important as the dynamics and quantum features, such as 
entanglement, of interacting many-body systems are further investigated. 
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